Class Web Page

R resource page


1. Effect of \(\bigtriangledown_s\) on Trend



1a. Linear Trend

  • Suppose \(s=12\), and your observation \(Y_t\) has linear trend, as well as the seasonality with period \(s=12\). \[ Y_t = a+b t + S_t + X_t. \] Then \[ \begin{align} \bigtriangledown_s Y_{t} &= \hspace3mm Y_{t+12}-Y_{t}\\\\ &= \hspace3mm a+b (t+12) + S_{t+12} + X_{t+12} \\\\ &- \hspace3mm a+b (t) \hspace15mm + S_t \hspace15mm + X_t \\\\ &= \hspace3mm b (12) + X_{t+12} - X_{t}\\ \end{align} \] This is statonary series with constant \(b(12)\). (called drift by auto.arima())


  • If your observation \(Y_t\) has quadratic trend, as well as the seasonality. \[ Y_t = a+b t + c t^2 + S_t + X_t. \] Then, \[ \begin{align} \bigtriangledown_s Y_{t} &= \hspace3mm Y_{t+12}-Y_{t} \\\\ &= \hspace3mm a+b (t+12) + c(t+12)^2 + S_{t+12} + X_{t+12} \\\\ &- \hspace3mm a+b (t) + c(t^2) \hspace5mm + S_t \hspace5mm + X_t \\\\ &= \hspace3mm b (12) + c(24 t) + c(12^2) + X_{t+12} - X_{t}\\ \end{align} \]


  • This is statonary series with linear trend \(b(12) + c(12^2) + c(24) t\).


  • If you take \(\bigtriangledown\) again, who will be left?



1b. Linear Trend and WN

  • Suppose \(s=12\), and your observation \(Y_t\) has linear trend, \[ Y_t = a+b t + S_t + e_t \hspace10mm \mbox{ where } e_t \sim N(0,1) \]

  • Taking \(\bigtriangledown_{12}\) yields, \[ \begin{align} Y_{t+12}-Y_{t} &= \hspace3mm a+b (t+12) + S_{t+12} + e_{t+12} \\\\ &- ( a+b (t) \hspace15mm + S_t \hspace5mm + e_t ) \\\\ &= \hspace3mm \hspace10mm b(12) + e_{t+12} - e_{t} \end{align} \]

  • What kind of process is this? \[ K_t \hspace3mm = \hspace3mm e_{t} - e_{t-12} \]

  • Linar regression with Seasonal average will be a better approach.



Example

## Series: diff(X, 12) 
## ARIMA(0,0,0) with non-zero mean 
## 
## Coefficients:
##         mean
##       6.1451
## s.e.  0.4902
## 
## sigma^2 estimated as 45.41:  log likelihood=-624.95
## AIC=1253.9   AICc=1253.96   BIC=1260.37

## Series: diff(X, 12) 
## ARIMA(0,0,0) with non-zero mean 
## 
## Coefficients:
##         mean
##       6.1451
## s.e.  0.4902
## 
## sigma^2 estimated as 45.41:  log likelihood=-624.95
## AIC=1253.9   AICc=1253.96   BIC=1260.37
## Series: diff(X2, 12) 
## ARIMA(0,0,0)(2,0,0)[12] with non-zero mean 
## 
## Coefficients:
##          sar1     sar2    mean
##       -0.6780  -0.2089  6.1699
## s.e.   0.0719   0.0747  0.2172
## 
## sigma^2 estimated as 29.7:  log likelihood=-586.82
## AIC=1181.63   AICc=1181.85   BIC=1194.58

##   B-L test H0: the sereis is uncorrelated
##   M-L test H0: the square of the sereis is uncorrelated
##   J-B test H0: the sereis came from Normal distribution
##   SD         : Standard Deviation of the series
##       BL15  BL20  BL25  ML15  ML20    JB   SD
## [1,] 0.363 0.152 0.115 0.107 0.173 0.608 5.42
## Series: diff(X2, 12) 
## ARIMA(0,0,0)(0,0,1)[12] with non-zero mean 
## 
## Coefficients:
##          sma1    mean
##       -1.0000  6.1078
## s.e.   0.0933  0.0683
## 
## sigma^2 estimated as 21.79:  log likelihood=-572.31
## AIC=1150.61   AICc=1150.74   BIC=1160.32

##   B-L test H0: the sereis is uncorrelated
##   M-L test H0: the square of the sereis is uncorrelated
##   J-B test H0: the sereis came from Normal distribution
##   SD         : Standard Deviation of the series
##       BL15  BL20  BL25  ML15  ML20    JB    SD
## [1,] 0.929 0.854 0.875 0.143 0.214 0.102 4.656



1c. Random Trend

  • Suppose \(s=12\), and your observation \(Y_t\) has random walk without drift as trend, together with the seasonality. \[ Y_t = W_t + S_t + X_t. \] Where \[ W_t \hspace3mm = \hspace3mm \sum_{i=1}^t e_i \hspace10mm e_i \sim_{iid} N(0,\sigma^2). \]

  • If you take seasonal difference, \[ \begin{align} \bigtriangledown_{12} Y_{t} &= \hspace3mm W_t \hspace5mm + S_t \hspace5mm + X_t \\ \\ &- \hspace3mm W_{t-12} - S_{t-12} - X_{t-12} \\ \\ &= \hspace3mm (W_t - W_{t-12}) + (X_t - X_{t-12}) \\ \end{align} \]

  • The first part, \[ W_{t} - W_{t-12} \hspace3mm = \hspace3mm \sum_{i=0}^{t-1} e_{t-i} - \sum_{i=12}^{t-1} e_{t-i} \\ \hspace3mm = \hspace3mm \sum_{i=0}^{11} e_{t-i} \hspace10mm \sim N(0,12\sigma^2) \]

  • This is MA(11) with unit root. (non-invertible)



Example

## Series: diff(RW2, 12) 
## ARIMA(1,0,0)(1,0,0)[12] with non-zero mean 
## 
## Coefficients:
##          ar1     sar1     mean
##       0.9458  -0.5853  -0.5674
## s.e.  0.0224   0.0590   0.9006
## 
## sigma^2 estimated as 1.301:  log likelihood=-293.38
## AIC=594.75   AICc=594.97   BIC=607.7

## Series: diff(RW2, 12) 
## ARIMA(0,0,12) with non-zero mean 
## 
## Coefficients:
##          ma1     ma2     ma3     ma4     ma5     ma6     ma7     ma8
##       0.9804  0.9841  1.0436  1.0251  1.0400  1.0394  1.0251  1.0419
## s.e.  0.0877  0.0927  0.1070  0.1184  0.1244  0.1202  0.1198  0.1100
##          ma9    ma10    ma11     ma12     mean
##       0.9816  0.9795  0.9610  -0.0382  -0.5658
## s.e.  0.1061  0.0964  0.1055   0.0826   0.7925
## 
## sigma^2 estimated as 0.9217:  log likelihood=-269.06
## AIC=566.11   AICc=568.54   BIC=611.42

## Series: diff(RW.d, 12) 
## ARIMA(1,0,0)(1,0,0)[12] with non-zero mean 
## 
## Coefficients:
##          ar1     sar1    mean
##       0.9458  -0.5853  5.4326
## s.e.  0.0224   0.0590  0.9006
## 
## sigma^2 estimated as 1.301:  log likelihood=-293.38
## AIC=594.75   AICc=594.97   BIC=607.7
## Series: diff(RW.d, 12) 
## ARIMA(0,0,11) with non-zero mean 
## 
## Coefficients:
##          ma1     ma2     ma3     ma4     ma5     ma6     ma7     ma8
##       1.0112  1.0141  1.0733  1.0555  1.0670  1.0670  1.0555  1.0733
## s.e.  0.0645  0.0685  0.0922  0.1024  0.1181  0.1094  0.1082  0.0910
##          ma9    ma10   ma11    mean
##       1.0141  1.0112  1.000  5.4364
## s.e.  0.0864  0.0689  0.072  0.8175
## 
## sigma^2 estimated as 0.9172:  log likelihood=-269.16
## AIC=564.33   AICc=566.42   BIC=606.4



2. Back to CO2 Data

Monthly CO2 levels at NW territories CANADA from Jan 1959 through Dec 1997.

## Series: co2 
## ARIMA(0,0,0)(0,1,1)[12] 
## 
## Coefficients:
##        sma1
##       0.375
## s.e.  0.072
## 
## sigma^2 estimated as 3.795:  log likelihood=-250.49
## AIC=504.98   AICc=505.08   BIC=510.56
## Series: co2 
## ARIMA(0,0,11)(0,1,0)[12] 
## 
## Coefficients:
##          ma1     ma2     ma3     ma4     ma5     ma6     ma7     ma8
##       0.7503  0.5927  0.6235  0.7815  0.8266  0.8043  0.8355  0.8013
## s.e.  0.1028  0.1221  0.1157  0.1170  0.1113  0.1368  0.1101  0.1189
##          ma9    ma10    ma11
##       0.6169  0.4684  0.7954
## s.e.  0.1208  0.1151  0.0955
## 
## sigma^2 estimated as 0.6845:  log likelihood=-150.54
## AIC=325.08   AICc=328   BIC=358.53
## Series: co2 
## ARIMA(0,0,11)(0,1,0)[12] with drift 
## 
## Coefficients:
##          ma1     ma2     ma3     ma4     ma5     ma6     ma7     ma8
##       0.3972  0.6207  0.2726  0.7195  0.4868  0.7640  0.4740  0.8061
## s.e.  0.1611  0.1247  0.1690  0.0886  0.1660  0.1241  0.1696  0.1164
##          ma9    ma10    ma11   drift
##       0.1842  0.5929  0.3300  0.1430
## s.e.  0.1883  0.1446  0.2036  0.0371
## 
## sigma^2 estimated as 0.6598:  log likelihood=-145.02
## AIC=316.04   AICc=319.48   BIC=352.28
## Series: co2 
## ARIMA(0,0,0)(0,1,1)[12] with drift 
## 
## Coefficients:
##          sma1   drift
##       -0.9998  0.1527
## s.e.   0.1653  0.0018
## 
## sigma^2 estimated as 0.6634:  log likelihood=-157.82
## AIC=321.64   AICc=321.85   BIC=330

  • Take \(\bigtriangledown_12\) and it has drift with sMA(1) close to 1.

  • This is a good sign that original co2 has \[ Y_t = W_t + S_t + e_t \]



3a. Linear Trend with WN

## 
## Call:
## lm(formula = Ds.co2 ~ time(Ds.co2))
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -2.22331 -0.69056 -0.01021  0.64085  2.36388 
## 
## Coefficients:
##                Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  -3.633e+03  5.112e+01  -71.07   <2e-16 ***
## time(Ds.co2)  1.817e+00  2.557e-02   71.07   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.9327 on 130 degrees of freedom
## Multiple R-squared:  0.9749, Adjusted R-squared:  0.9747 
## F-statistic:  5052 on 1 and 130 DF,  p-value: < 2.2e-16

## Series: Reg1$resid 
## ARIMA(1,0,0) with zero mean 
## 
## Coefficients:
##          ar1
##       0.5273
## s.e.  0.0739
## 
## sigma^2 estimated as 0.6228:  log likelihood=-155.71
## AIC=315.42   AICc=315.51   BIC=321.19

##   B-L test H0: the sereis is uncorrelated
##   M-L test H0: the square of the sereis is uncorrelated
##   J-B test H0: the sereis came from Normal distribution
##   SD         : Standard Deviation of the series
##       BL15 BL20  BL25  ML15  ML20    JB    SD
## [1,] 0.386 0.13 0.074 0.691 0.817 0.416 0.789
## Series: Resid.co2 
## ARIMA(2,0,0)(2,0,0)[12] with zero mean 
## 
## Coefficients:
##          ar1     ar2    sar1    sar2
##       0.3858  0.2620  0.2460  0.2002
## s.e.  0.0852  0.0891  0.0901  0.0924
## 
## sigma^2 estimated as 0.5589:  log likelihood=-148.2
## AIC=306.4   AICc=306.88   BIC=320.82

##   B-L test H0: the sereis is uncorrelated
##   M-L test H0: the square of the sereis is uncorrelated
##   J-B test H0: the sereis came from Normal distribution
##   SD         : Standard Deviation of the series
##       BL15  BL20  BL25  ML15  ML20    JB    SD
## [1,] 0.758 0.425 0.487 0.641 0.602 0.904 0.739
## Series: Resid.co2 
## ARIMA(2,0,0)(1,0,0)[12] with zero mean 
## 
## Coefficients:
##          ar1     ar2    sar1
##       0.3970  0.2333  0.2862
## s.e.  0.0867  0.0916  0.0952
## 
## sigma^2 estimated as 0.5791:  log likelihood=-150.46
## AIC=308.92   AICc=309.24   BIC=320.45

##   B-L test H0: the sereis is uncorrelated
##   M-L test H0: the square of the sereis is uncorrelated
##   J-B test H0: the sereis came from Normal distribution
##   SD         : Standard Deviation of the series
##       BL15  BL20  BL25  ML15  ML20    JB    SD
## [1,] 0.547 0.277 0.238 0.373 0.408 0.847 0.755



Model for CO2 data

  • \(s=12\) because freq=12 was set in ts object.

  • Took \(\bigtriangledown_{12}\) then checked for sMA(1) paramter and MA(11) parameters.

  • The drift term (mean after the seasonal difference) was significant. With the drift, sMA(1) parameter was too close to 1, indicating that original series has linear trend with WN \(e_t\).

  • With observation \(Y_t\), \[ Y_t \hspace3mm = \hspace3mm L_t + S_t + X_t \] where \(X_t\) was \(ARMA(2,0) \times (1,0)\) \[ X_t = \phi_1 X_{t-1} + \phi_2 X_{t-2} + \Phi_1 X_{t-12} + e_t \]




3. Forecasting sARMA


  • Suppose we have ARIMA\((1,0,1) \times (1,1,0)_{12}\) model, \[ (1-\phi_1 B) (1-\Phi_1 B^{12}) \bigtriangledown_{12} Y_t = (1-\theta_1 B) e_t \]

  • This is same as modeling \(\bigtriangledown_{12} Y_t\) with ARMA(13,1), \[ X_t \hspace3mm = \hspace3mm \bigtriangledown_{12} Y_t \\ \Big(1-\phi_1 B - \Phi_1 B^{12} + \phi_1 \Phi_1 B^{13} \Big) X_t \hspace3mm = \hspace3mm (1-\theta_1 B) e_t \]

  • We know how to get a prediction \(\hat X(h)\) for ARMA(\(p,q\)), Our predictor for \(Y_t\) will be \[ \begin{align} Y_{n+1} &= \hspace3mm Y_{n+1-12} + \hat X(1) \\ Y_{n+2} &= \hspace3mm Y_{n+1-12} + \hat X(2) \\ &\vdots \end{align} \]




4. Tests for Seasonality


How do you check if you need to take a seasonal difference, instead of regular difference?

  • OCSB test (\(H_0\): Seasonal Unit Root Exists)

    • Default method in auto.arima()
    • Osborn-Chui-Smith-Birchenhall (1988) \
  • CH test (\(H_0\): Deterministic Seasonality Exists)

    • Canova-Hansen (1995)
  • Fore more detail: http://robjhyndman.com/hyndsight/forecast3/



4a. OCSB test

Osborn-Chui-Smith-Birchenhall (1988) test for

  • Default choice in for choosing value of \(D\).
    \[ \left\{ \begin{array}{ll} H_0: \mbox{ Seasonal Unit Root Exists}\\ H_A: H_0 \mbox{ is false} \\ \end{array} \right. \]

  • Small p-value means “Take Seasonal Difference”.





5. Summary

  • If you take \(\bigtriangledown_{12}\) of Linear trend + Seasonality, it eliminates the seasonality, and leaves constant \(12(slope)\).


  • But if you do the same with (Linear trend + Sesonality + WN), then it turns into \(12(slope) + (non-invertible sMA(1))\). Which is troublesome. Force fit sMA(1) to find out.


  • Seasonality and trend should be taken care of via regression and seasonal average, not by \(\bigtriangledown_{12}\).


  • If you take \(\bigtriangledown_{12}\) of Random Trend + Seasonality, then it leaves 12(drift) + (non-invertible MA(11)). Force fit MA(11) to find out.



6. R-code Only